import numpy as np
import Mixing_rules as MR
from Pure_VLE_TB import T_BOIL_CAL
from EOS_collection import CUBIC_EOS
from scipy.optimize import minimize

R = 8.314

#输入对应态参数
Tc = np.array([[364.8, 408.1]])
Pc = np.array([[4.610e6, 3.648e6]])
w = np.array([[0.148, 0.176]])
Z = np.array([[0.275, 0.281]])
M = np.array([[44, 44]])
Tb = np.array([[322.02, 373.42]])

#输入要计算的状态参数

x1 = np.array([[0.2]])
x2 = 1 - x1
x = np.concatenate((x1, x2))
y = np.zeros_like(x)
T0 = np.ones_like(x) * 300
T = np.sum(x.flatten() * Tb.flatten())
P = 2026.5e3

EOS = CUBIC_EOS(Pc, Tc, Z, w)

def BUBBLE(T, y):

    # 初始化
    ynew = x.flatten() * Pc.flatten() / ( np.sum ( x.flatten() * Pc.flatten() ) )
    ynew = ynew.reshape(-1, 1)

    # 计算状态方程参数
    a = EOS.PR_a(T)
    b = EOS.PR_b()

    # 计算液相逸度系数
    aL = MR.a_mix(a, x)
    bL = MR.b_mix(b, x)
    VL = 2 * bL
    VLnew = EOS.PR_V(P, VL, T, aL, bL)
    while np.sum ( np.abs(VLnew - VL) ) / np.sum ( VL ) > 1e-5:
        VL = VLnew
        VLnew = EOS.PR_V(P, VL, T, aL, bL)
    ZML = P * VL / R / T
    xa = MR.xya_mix(a, x)
    phiL = EOS.PR_PHI_MIX(P, VL, T, aL, bL, xa, b, ZML)

    # 求气相的组成
    while np.sum(np.abs(ynew - y)) > 1e-5:
        # 计算气相逸度系数
        y = ynew
        aV = MR.a_mix(a, y)
        bV = MR.b_mix(b, y)
        VV = R * T / P
        VVnew = EOS.PR_V(P, VV, T, aV, bV)
        while np.sum ( np.abs(VVnew - VV) ) / np.sum ( VV ) > 1e-5:
            VV = VVnew
            VVnew = EOS.PR_V(P, VV, T, aV, bV)
        ZMV = P * VV / R / T
        ya = MR.xya_mix(a, y)
        phiV = EOS.PR_PHI_MIX(P, VV, T, aV, bV, ya, b, ZMV)
        ynew = x.flatten() * phiL / phiV
        ynew = ynew.reshape(-1, 1)
    
    return (np.sum(y) - 1) ** 2

T = minimize(BUBBLE, T, y, method = 'L-BFGS-B').x
print(f"泡点温度是{T}K")